#include <istream>
#include <vector>
#include <string>
#include <fstream>
#include <iostream>
#include <algorithm>
#include <map>
#include <pthread.h>
#include <boost/format.hpp>
#include <boost/lexical_cast.hpp>
#include <boost/filesystem/path.hpp>
#include <cstdlib>

#include <print_error.h>
#include <tree_UPGMA.h>
#include <tree.h>
#include <tree_ops.h>
#include <portable_timer.h>
#include <progress_indicator.h>
#include <sparse_matrix.h>
#include <sparse_matrix_io.h>
#include "get_arg.h"
#define DEBUG_LEVEL 6
#include <debug_func.h>

using boost::format;
using boost::lexical_cast;
using boost::filesystem::path;
using std::vector;
using std::ostringstream;

using std::string;
using std::fill;
using std::ofstream;
using std::cout;
using std::cin;
using std::ifstream;
using std::cerr;
using std::map;
using std::pair;
using std::istream;
using std::max;
using std::min;
using namespace bcpp_matrix;

struct t_incorrect_pairs
{
	string cluster_name;
	unsigned cnt_genes;
	unsigned cnt_missing_pairs;
	unsigned cnt_expected_pairs;
	t_incorrect_pairs(const t_sparse_matrix_f& data, const string& cluster_name_):
		cluster_name(cluster_name_),
		cnt_genes(data.size()),
		cnt_missing_pairs(data.get_count_missing()),
		cnt_expected_pairs(data.get_count_pairs())
	{}
};


//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	main

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

int main (int argc, char *argv[])
{
	// we use cout to print progress so flush per call
//	cout.setf(std::ios_base::unitbuf);

	std::ios::sync_with_stdio(false);

	t_program_args args;
	try
	{
	if (!process_command_line_options(argc, argv, args))
		return 1;




	unsigned cnt_cluster = 0;
	unsigned cnt_subclusters = 0;
	//
	//	save clusters with missing entries
	//
	vector<t_incorrect_pairs> problems;



	//
	//	Multiple clusters can be read in with pairwise input format
	//
	args.VERBOSE(1) << "\tBuilding UPGMA trees:\n";
	while (cin)
	{
		++cnt_cluster;


		//
		//	read data from input file
		//
		vector<t_sparse_matrix_f>	data;
		args.VERBOSE(2) << "\tReading pairs and single linkage clustering...\n";
		{
			t_sparse_matrix_f initial_matrix;

			t_progress_indicator dots(args.VERBOSE(1), 20000);
			t_auto_progress_indicator auto_dots(dots, args.is_verbose() == 2);
			dots.set_every(200000);
			bcpp_matrix::read_matrix(cin,  initial_matrix,  args.input_format, dots);
			if (initial_matrix.matrix_name.empty())
				initial_matrix.matrix_name = "Cluster# " + lexical_cast<string>(cnt_cluster);
			dots.print_indicator(':');
			initial_matrix.get_single_linkage_clusters(data);

		}
		if (!data.size())
			break;

		// get how many things we have to resolve for time estimate
		unsigned cnt_matrix_cells_to_cluster = 0;
		for (unsigned i = 0; i < data.size(); ++i)
		{
			if (data[i].get_count_filled() == 0)
				continue;
			cnt_matrix_cells_to_cluster += data[i].size() * data[i].size() * 2;
		}


		t_progress_indicator dots(args.VERBOSE(1), 100000);
		dots.set_limit(dots.count_dots + cnt_matrix_cells_to_cluster);
		dots.set_every(std::max(100000U, (dots.count_dots + cnt_matrix_cells_to_cluster) / 500));
		dots.use_timer();
		for (unsigned i = 0; i < data.size(); ++i)
		{
			if (data[i].get_count_filled() == 0)
				continue;

				++cnt_subclusters;

			// if has multiple sub-clusters, append index of subcluster after period
			if (data.size() != 1)
				data[i].matrix_name += " Subcluster# " + lexical_cast<string>(i);

			t_sparse_matrix_f& sub_matrix = data[i];
			// save problems
			if (sub_matrix.get_count_missing())
				problems.push_back(t_incorrect_pairs(sub_matrix, data[i].matrix_name));

			string header = ">" + sub_matrix.matrix_name;
			args.VERBOSE(2) << "\t" << header << "\n";
			//args.VERBOSE(2) << "\t" << sub_matrix.size() << " genes\n";


			bcpp_tree::t_tree upgma_tree;

			//888888888888888888888888888888888888888888888888888888888888888888888888

			// UPGMA

			//888888888888888888888888888888888888888888888888888888888888888888888888
			{
				args.VERBOSE(2) << "\tConstruct UPGMA tree...\n";
				t_auto_progress_indicator auto_dots(dots, args.is_verbose() == 2);
				upgma_tree::upgma_make_tree(upgma_tree, sub_matrix,
											false /*args.fill_lacunae*/, dots);

			}


			// print new_hampshire tree
			{
				if (args.header)
					cout << header << "\n";
				upgma_tree.assert_consistency();
				cout << upgma_tree << std::endl;
				if (args.add_matrix_to_output)
				{
					if (!args.header)
						sub_matrix.matrix_name = "";
					dots.set_every(100);
					t_auto_progress_indicator auto_dots(dots, args.is_verbose() == 2);
					cout << t_sparse_matrix_ostrm_wrapper<float>(sub_matrix, args.output_format,  
														  dots);
				}


//#if (!defined(NDEBUG))
//				bcpp_tree_ops::do_calculate_branch_lengths branch_lengths(upgma_tree);
//				//std::vector<eNODETYPE>	node_types(tree.get_node_count(), eNONE);
//				bcpp_tree_ops::op_print_tree_elaborate
//											print_tree_elaborate
//											(upgma_tree,
//											args.f_tree,
//											branch_lengths.lengths_to_leaves,
//											branch_lengths.lengths_to_root);
//				print_tree_elaborate.print();
//#endif
			}

		}
		if (args.input_format != bcpp_matrix::eMIO_PAIRS)
			break;
	}

	//
	//	print clusters with missing entries
	//
	if (problems.size())
	{
		std_print_warning(args.VERBOSE(),  "Some matrixes had missing pairs.");
		for (unsigned i = 0; i < problems.size(); ++i)
		{

			args.CERR() << format("\t\t%|1$4|: "
									"%|2$4| genes, "
									"%|3$4| gene pairs expected, "
									"%|4$4| gene pairs missing.\n ")
					% problems[i].cluster_name
					% problems[i].cnt_genes
					% problems[i].cnt_expected_pairs
					% problems[i].cnt_missing_pairs;
		}
	}
	args.VERBOSE() << "\tAll finished!!\n\n";
	}
	catch(std::exception& e)
	{
		std_print_error(args.CERR(), e.what());
		return 1;
	}

	return 0;
}

